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ABSTRACT 

The behavior of a fully-submerged hydrofoil has been 
studied by using Root Locus technique. A multivariable con- 
cept was used in designing the compensators of the system. 

A digital computer simulation program using Digital Simula- 
tion Language or DSL is produced to study the performance on 
the digital computer IBM system 360/67. The Highpoint 
PC(H)-1 is used as the model of this simulation. 
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I. INTRODUCTION 



A. FOIL TYPES 

The main idea of the hydrofoil craft is to reduce the 
resistance of the craft by reducing the area in contact with 
the water and thus increase the speed for the same power. 

One way to accomplish this is fitting foils below the hull 
of a craft, so that the craft is lifted clear of the water 
at high speed. There are two principal foil arrangements in 
use today, classed as surface piercing foil and fully-sub- 
merged foil, as shown in Figure 1. 



The surface-piercing foil is supported by vertical 
struts and the lifting surface extends in a V-shape from the 
water surface down through the water to the point of conver- 
gence of the V. The surface-piercing foil provides a 




Surface-Piercing Foil 



Fully- Submerged Foil 



FIGURE 1 



TYPE OF HYDROFOIL CRAFT 
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self-stabilization for the craft in pitch, height and roll. 
The self-stabilizing aspect of the surface piercing foil is 
easy to see. Suppose the craft pitches up slightly, then 
less area of the forward foil and more area of the after foil 
is submerged. This results in less lift from the forward 
foil and more lift from the after foil so the craft will 
settle back again, the same stabilization occurs in roll. 

For the height as the craft speed increases the whole craft 
rises up until the amount of wetted surface of the foils is 
just sufficient to create a lift equal to the weight of the 
craft . 

The fully-submerged foil cannot adjust itself to chang- 
ing lift demands by a change in area, but only by a change 
in incidence, either of the whole foil or of attached flaps. 
This involves the use of some sensing arrangement to control 
the incidence by monitoring the height of the craft above 
water surface. The fully-submerged foil is a more sophisti- 
cated design, and the control mechanism is likely to be com- 
plicated. But the fully-submerged foil has a very distinct 
advantage when operating in heavy seas. The surface-piercing 
foil will respond to even relatively small changes in the 
surface of the sea and give a very rough ride. The fully- 
submerged foil, on the other hand, is much less affected by 
the change in the surface of the sea. This paper will deal 
only with the fully-submerged foil craft. 
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B. FOIL ARRANGEMENT 



Most of the fully-submerged foil craft have two foils, 
one at the forward and the other at the aft. Both foils are 
supported by the vertical struts. Figure 2 shows the con- 
figuration of main foil and secondary foil; Figures 3A and 
3B show foils in canard arrangement and in airplane arrange- 
ment respectively. 

The simulation in this paper is based on the canard ar- 
rangement, which has a secondary foil as a forward foil and 
a main foil as the aft foil. The HIGH POINT PC(H)-1 was the 
model for this simulation. Drawing of the craft and the 
foils are shown in Figures 4A and 4B. 
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FIGURE 2 

FOILS CONFIGURATION 
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FIGURE 3A 

FOIL IN CANARD ARRANGEMENT 




FIGURE 3B 

FOIL IN AIRPLANE ARRANGEMENT 
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\^>.(SYMETRICAL) 

FIGURE 4A 

DRAWING OF FOIL AND STRUT CONFIGURATION 



Aft Foil Fq/waxp 
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DRAWING CF FOIL AND FLAP CONFIGURATION 



II. HYDROFOIL EQUATIONS STUDY 



A. COORDINATE SYSTEMS 

Three coordinate systems will be used in this study of 
hydrofoil craft dynamics. They are: the Body Axes, the 
Water Axes and the Earth Axes. Each is in a right hand 
orthogonal system. 

The body axes coordinate system (x,y,z) shown in Figure 

* 

5 has its origin at the center of gravity of the craft and 
is fixed relative to the craft. This system is initially 
coincident with the origin of the earth axes system, however, 
it moves with the craft as time progresses. X is positive 
forward, Y is positive to starboard and Z is positive down- 
ward . 




FIGURE 5 



ORIENTATION OF THE BODY AXES WITH RESPECT TO THE 



CRAFT AND DIRECTIONS OF POSITIVE VELOCITIES 
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The earth axes coordinate system (Xg,Yg,Zg) is fixed 
relative to the earth surface. The origin may be chosen at 
any convenient location, provided the body system origin is 
initially at the same point. The X^, axis lies in the hori- 
zontal plane; it is initially coincident with, and positive 
in the same direction as the body X axis. The Y^ axis lies 
in the horizontal plane and is positive to starboard when 
the observer is facing the positive X^ direction. The 
axis is normal to the horizontal and is positive downward. 
Water surface motion and craft motion relative to calm water 
are described in this coordinate system. 

The water axes coordinate system is aligned with the 
relative velocity vector and resolves into "Lift," "Drag," 
and Side -Force" direction. 




It will be specified that the order of rotation from 
water axes to body axes is: 

1) The (D,S) axes rotate about the L-axis through the 
angle $ to the position (D^.S^) where is coincident with 
the body Y-axis. 

2) The (D^,L) axes rotate about the S^-axis, or more 
pertinently about the Y-axis, through the angle a to the 
position (Dg.L^) which corresponds to the (-X,-Z) axes. 



B. EQUATION OF MOTION 

The development of the six equations of motion for a 
hydrofoil follows the same procedure as those for a displace- 
ment type hull. This development is clearly outlined in 
Ref. 2 and 8 and therefore will not be done here. 

The equations of motion are: 

Translation 



u = -F„ - QW + RV 
m X 

v = -F v + PW - RU 
m Y 

w = ^F + QU - PV 
m z 

Rotation 

p = I^ [L - QR < I zz- 1 Yy> + ( f!+C!P > I xz 1 

Q = i^[ m - rp ( I xx" I ZZ )-(P2_R2)I XZ ] 

r = I i-[N-PQ(I Yy -I xx )-(QR-P)I xz ] 



( 2 - 1 ) 

(2-2) 

(2-3) 

(2-4) 

(2-5) 

(2-6) 



The terms RV, PV, and PW are very small and may be neg- 
lected. Similarly, the product of inertia quantity, lyz’ i s 
less than 10% of Iy X and less than 4% of I^Z' Consequently, 
when considered in conjunction with the magnitudes of (R+PQ), 



(P z -R 2 ) and (QR-P), the terms containing I vz may be neglected 



for this craft. 

Equations (2-1) through (2-6) are linear and rotational 
accelerations in body axes. The velocity terms in body axes 
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can be obtained directly by integrating each equation one 
time. 



u = 


/ 


Udt 


(2-7) 


V = 


/ 


Vdt 


(2-8) 


w = 


/ 


Wdt 


(2-9) 


p = 


/ 


Pdt 


(2-10) 


Q = 


/ 


Qdt 


(2-11) 


R = 


/ 


Rdt 


(2-12) 



Linear displacements and angular positions with respect 
to body axes could be obtained by an additional integration, 
but such quantities are of limited usefulness. The dis- 
placements and angle usually of interest are those with 
respect to earth axes. In order to calculate the craft's 
position in earth axes coordinates, the craft velocity vec- 
tor must be transformed from body axes coordinates into 
earth axes coordinates as follow: 

Ug = Ucos6cos^+V[cosi|Jsin0siniJj-sini|>cos(J>] 

+W[cos\Jjsin0cos<l>+siniJjsin<j> ] (2-13 ) 

V E = Ucos0simJj+V[cosi|JCos<}>+siniJ;sin0sin<j>] 

+W[ sinijjsin0cos(j>-cosijjsin<j>] (2-14 ) 

= -Usin0+Vcos0sin<f>+WcosOcos<}> (2-15) 

Only the craft's velocity in the downward, Z^., direction 
is required to calculate its height above the water surface. 
That is, the vertical position in earth axes follows from 

Z £ = / W E dt ( 2-1G ) 
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The craft's position coordinates are obtained from 

X E = / U E dt (2-17) 

Y E = / V £ dt (2-18) 



The Euler angles, 4> , 0, and ip are also of interest and 

• • • 

can be found by integration of Euler rates <j>, 0, and ip. The 
rates in turn must be derived from the body angular rates. 
The Euler rates are not easily obtained because they occur 
about axes which are not orthogonal. This fact can be ap- 
preciated by referring to Figure 6 and recalling how the 
angles were defined, ip is a rotation about the Z 0 axis, 0 
is a rotation about the Y x axis, and <J> is a rotation about 
the X axis. They result in the following set of equations: 



<j> = P+tjjsin0 


f o -» r\ \ 


• 

0 = Qcos<j>-Rsin<j> 


(2-10) 


• © 

ip = (Qsin<J)+Rcos<}> )cos0 + (<j>-P)sin0 


(2-21) 



Studies of this craft show that even under the most 
drastic abnormal situations, pitch angles are not expected 
to reach 10° and under normal operating conditions, roll 
angles will not exceed 3°. Therefore, the following assump- 
tions are justified; 



sin0 = 0 
COS0 = 1 



( 2 - 22 ) 



The Euler angles are obtained by integration of equations 
(2-19) through (2-21). 
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FIGURE 6 

DEFINITION OF POSITIVE DIRECTIONS OF EULER ANGLE 



R6 



0. CALCULATION OF VELOCITY COMPONENTS 

The total velocity components of each foil in body axes 
can be expressed in terms of craft linear and angular veloc 
ities by considering the craft geometry as represented in 
Figure 7 and Table I . 

The resultant equations for the example shown are: 



Center Foil 

U C = u + l zcf' q 

V C = v " l zcf* p + l xcf’ r 

w c = w - l xcf *Q 

Port Foil 

• Up = u + l zpf q ' l ypf - r 

Vp = v - Lgpj- P + l xpf' R 

Wp - W + L YPF - p - l xpf - q 



(2-23) 

(2-24) 

(2-25) 

(2-26) 

(2-27) 

(2-28) 



Starboard Foil 



S 

J S 

!1 s 



= u + l Z sf* q " l ysf* r 


(2-29) 


= v - L ZSF* p + l xsf* r 


(2-30) 


= W + Ly SF -P - L XSF* Q 


(2-31) 



Mid-Foil Segment 



U M = U + L ZMF* Q 

V M = v - L ZMF' P + L XMF 



R 



W M = W - L 



XMF 



Q 



(2-32) 

(2-33) 

(2-34) 
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FIGURE 7 

DEFINITION OF SYMBOLOGY 

Point of application of foil forces 
Point of application of strut forces 



l xcf 


L xcs 




l xpf 


L XPS 


l xsf = L xss = l xmf 


^SF 


hfSS 


l ypf = SrPS 


l zpf 


l zsf 


l zmf 


L ZPS 


and L zss 


are variables and in 
equal . 
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Point of Symbol for Force Along Coordinate From CG to Point 

Application of or Parallel to of Application Relative to 
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of Prop. Thrust 



These foil velocities are in body axis terms and can be 
used directly to calculate the angles of attack and angles 
of side slip. 

D. ANGLE OF ATTACK AND ANGLE OF SIDE SLIP 

The angle of attack of the foil is shown in Figure 8A 
and defined as the angle whose tangent is the total relative 
velocity between the water and the foil in the body Z direc- 
tion, divided by the total relative velocity between the 
water and foil in the body X direction 

W 

a = arc tan ■— (2-35) 

u r 

If the foil is not aligned with the body X-axis, then 
the foil- has a fixed angle of attack and the total angle of 
attack becomes 

a TOTAL = a FIXED + a (2-36) 

The side slip angle is shown in Figure SB and can be 
represented by 

V ri 

$. = arc tan ' (2-37) 

1 Aj 2 . +v 2 . +w 2 . 

n ri ri 

During actual craft operations, the angles will be less 
than 10°, and will permit use of the small angle approxima- 
tion. Also the V and V. 7 velocities can be considered very 
small in relation to the U velocity. Applying these simpli- 
fications, the equations for a and 3 become: 
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• FIGURE 8a 
ANGLE OF ATTACK 
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FIGURE 8b 

ANGLE OF SIDE SLIP 
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E. EXPANSION OF FORCE AND MOMENT TERMS 

With the aid of Figure 7 and Table I, the forces and 
moments are set down. It is important that when a numerical 
value is inserted for one of the alphabetical dimensions, 
the sign of the value must be positive or negative according 
to positive and negative directions from the axis system 
origin at the C.G. Also note that horizontal foils make no 
contribution to the body Y force and vertical struts do not 
contribute to the body Z force. One of the assumptions made 
upon commencing simulation was that the struts contributed 
zero lift and that the foils contributed zero side force. 
Therefore, all body Y motion originates at the struts and 
all body Z motion originates at the foils. 

The force and moment equations become: 



F = F +F +F +F +F +m +T 
X XCF XPF r XSF XCS XSS gx X 



+m +T 



( 2 - 38 ) 




( 2 - 39 ) 



F = F +F +F +F +m 
Z ZCF ZPF ZSF ZMF gz 



( 2 - 40 ) 



L = (F 




( 2 - 41 ) 






( 2 - 42 ) 



N = (F v L v ) r .„+(F v L Y )_+(F v L v ) 



Y X'CS v Y X 'PS v Y X'SS 



( 2 - 43 ) 
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All of these forces have a hydrodynamic origin except 
for the thrust and gravity terms. The thrust is always 
associated with the body axis and requires no further expan- 
sion. The gravity terms can be simplified using small angle 
approximation for the Euler angles (see equations A-31,32,33). 

The general expression for hydrodynamic forces in water 
axes coordinates is 

F. = qACL i = L , D , S (2-44) 

where L,D,S represent lift, drag and side force respectively. 
The lift, drag, and side force coefficients very as func- 
tions of angle of attack, angle of side slip, submergence, 
velocity, flap elevator and rudder positions. 

Ideally, mathematical expressions would have been de- 
veloped to correctly depict the interrelation of all the 
variables which affect the force coefficients. However, 
hydrofoil technology has not advanced to a position which 

would yield such expressions. Some of the major problems 

* 

encountered in deriving a mathematical expression are: 

1) The occurrence of cavitation and ventilation which 
are completely unpredictable. 

2) Even more unpredictable is the cessation of cavita- 
tion and ventilation. 

3) The nonlinearity of the hydrodynamic coefficients. 

4) Lack of sufficient test data to completely describe 
all of the above. 
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F. SPECIFICATION OF UNITS 



Considerable confusion resulted from the interpretation 
of the various references concerning the units used in de- 
fining Euler angles and angular rotations. This problem has 
been alleviated here by the specific definitions in the 
Table of Symbols. 
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III. HYDROFOIL SIMULATION 



A. SIMULATION OBJECTIVES 

The object of the simulation was to commence the runs 
with the craft foilborne and in a steady state condition. 
Once this condition was achieved, the step response of the 
craft would be studied by applying step functions to excite 
motion along and about each body axis. That is, perturba- 
tions should be introduced to separately excite: 

1) Motion along the X-axis 

2) Motion about the X-axis 

3) . Motion along the Y-axis 

4) Motion about the Y-axis 

5) Motion along the Z-axis 

6) Motion about the Z-axis 

After the step response was obtained, a sinusoidal sea 
would be inserted to observe the craft's motion in a seaway 
The following assumptions were made prior to beginning simu 
lat ion : 

1) Craft equations of motion are valid only for the 
foilborne mode. 

2) Weight of the craft remains constnat. 

3) Hydrodynamic coefficients are based on fully wetted 
surfaces, i.e. no cavitation or ventilation. 

4) No constraint was placed on craft's heading. 
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5) Foils are considered flat surfaces vice dihedral/ 
anhedral. 

6) The craft may be operated in calm water or seaway. 

B. COMPUTER SIMULATION 

The computer language used in this simulation was Digital 
Simulation Language or DSL, and based on the computer pro- 
gram in Ref. 1. 

Before the main program could be assembled, subroutines 
had to be written to find the hydrodynamic coefficients of 
foil lift and drag, and strut drag and side slip. Using the 
simulation curves in Ref. 3 as a guide, the curves were 
stored as data points in the setup routine. The subroutines 
"INTERP" and "INTERP1" are interpolation routines to obtain 
the proper coefficient from the curves. For INTERP, values 
of angle of attack and flap angle are used to enter the 
curves and obtain a value of foil drag. Angle of side slip 
and submergence are used to enter INTERP1 to obtain strut 
drag. The curve for the lift coefficient of a foil in the 
fully wetted region is straight line so this coefficient is 
found by merely solving the equation for the straight line. 
Curves plotted from the subroutines are shown in Appendix B. 
Listings of the subroutines are shown in the computer pro- 
gram section. The main program was then assembled using 
Figure 9 as a reference for data flow. One part before the 
end of the main program was opened for an automatic control 
block. A complete listing of the main program is shown in 
the computer program section. 
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BLOCK DIAGRAM OF EQUATION RELATIONSHIPS 




C. AUTOMATIC CONTROL 



Control signals originate at the pilot's controls, at the 
motion sensors, and at the position transducer. These are 
combined and processed in the electronics to produce the 
signal for the servo- valve. The pilot will be able to move 
the helm, a lever, or a knob to introduce steering, altitude 
and attitude commands, the motion sensors sense errors be- 
tween the commanded and actual craft motions and produce 
electrical signals proportional to these errors. The posi- 
tion transducer provides the feedback signal so that the 
control surface stope moving when it has reached the position 
corresponding to the processed error signals. 

It is necessary to consider how the control system com- 
ponents affect the equations of motion through the control 
surface deflections. In general, a signal originating at 
the pilot or at the motion sensors will be modified by the 
individual dynamics of the components in the signal path be- 
fore the control surface motion actually occurs. Motion 
sensors to be discussed are vertical gyros, rate gyros, ac- 
celeromaters and height sensors. 

1. Vertical Gyros 

A vertical gyro of the type of interest here will 
produce one signal directly proportional to the pitch angle 
of the craft and one signal directly proportional to the 
roll angle of the craft. 
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2. Rate Gyros 



A rate gyro will sense an angular rate about one of 
the three body axes. (Three-gyro "rate packages" are avail- 
able if required for a particular application.) The dynam- 
ics of the gyro are such as to produce an output signal 
which is related to the actual angular rate by a second order 
differential equation. Using the roll rate P, as an example, 
the equation is 



P = 



0) 



2 W 

n 



2£o> ~ + U) 2 P) 

n dt n 



(3-1) 



or in transfer function form for all rotational rates: 

P = Q = R = ^ 

P Q R S 2 +2<> S+o) 2 

n n 

where 

P,Q,R = Outputs of rate gyros (RAD/SEC) 

P,Q,R = Craft rotational rates (RAD/SEC) 
a) = Natural frequency of a particular gyro 

(RAD /SEC) 



C = Damping ratio of a particular gyro 

(Dimensionless) 

S = Laplace transform variable 

3 . Accelerometer 

Outputs of accelerometers are also assumed to relate 
to the craft motions by second order differential equations. 
In general form, the equation is: 



^ d 2 a . aa . 

< ” + 2cu n Jjr 1 + u’a. ) 



da . 



a 



o)j 2 dt 2 



n i 



(3-3) 
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or : 




a 



i 



0) 



n 



S 2 +2Cw S+w 2 
n n 



(3-4) 



where i will become a symbol corresponding to the location 
of a particular accelerometer. 

To identify the accelerations sensed by the instru- 
ments, it will be assumed that there is one accelerometer at 
the C.G. mounted such as to measure lateral acceleration, 
and one accelerometer above each foil mounted in the X-Y 
plane so as to measure accelerations parallel to the Z-axis. 

The lateral acceleration at the C.G. is: 

F 

o = y ( 



Taking angular acceleration into account , the accel- 
erations in the Z-direction at points above a center, port, 
and starboard foil are: 



ZC 



F 

_Z T n 

m L XCF Q 



A = — 
A ZP 



m + L XPF* Q L YPF* P 



A ZS m + ^SF^ + S^SF*^ 



(3-6) 

(3-7) 

(3-8) 



4. Height Sensor 

Ultrasonic and sonar types are either in use or ex- 
pected to be used. For this study it will be assumed that 
the dynamic characteristics of the sensor, whatever the 
type, will be such as to give a direct indication of the in- 
stantaneous height of the sensor above the water. 
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IV. TRANSFER FUNCTION AND ROOT LOCUS DESIGN 



A. PITCH-HEAVE MOTION TRANSFER FUNCTION 

The best way to study the stability of the craft was to 
find the transfer functions of the craft and use either Root 
Locus plots or Bode plots to be an aid in a stability study 
and also in designing the compensators of the system. 

The transfer functions of the craft can be derived di- 
rectly from the equations of motion. In the pitch-heave 
mode, equations (2-3) and (2-5) were the equations that de- 
scribed the motion in the Z-axis and motion about the Y-axis. 
After neglecting the small terms, equations (2-3) and (2-5) 
were as follow 



YY 

From equations (2-35) and (2-36) the angle of attack of 
each foil can be written as follows: 




(3-9) 



Q = _L_ -M 



(3-10) 



W 




a 



C 



(3-11) 



W 




a 



S 



(3-12) 



W 




a 



P 



(3-13) 



W 




a 



M 



(3-14) 
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W RP W + ^PF* P L XPF* Q 



(3-15) 



\s ■ 


= w + w p - w Q < 3 - i6 > 


W RM 1 


= W + ^lF -9 (3 ' 17) 


ff RC ' 


■ W - L xcf -Q (3-18) 



Assume there is no significant coupling from the roll, 
yaw, sway or surge motions. Equations (3-15), (3-16), (3-17) 
and (3-18) will be 



w - 

RP 


= W - L XPF -Q 


W 

W RS 


■ * - Lxsf^ 


W RM 


II 

1 

t-* 

<o 


W RC 


■ * - l xcf* ^ 



And equations (3-11), (3-12), (3-13) and (3-14) will be 



a c 


w l xcf at 

u u al f 


a s = 


w l xsf n , , T 
U u + AL s 


p 

s 

!! 


w l xmf n AT 
U D ' Q + AL 11 


p 

II 


w l xpf _ 
u U ^ AL P 



^SF L XPF L XMF 
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Let ALp , ALg ALp AL^j ^2 



l xcf . 

u 



Q - Qi 



^SF 

U 



Q — Q2 




a c = W -Q1+A1 

a S = W '-Q 2 +A 2 " a p * a M 



From Equation (2 40) F z F zCF +F ZPF +F ZSF +F ZMF +F gZ 

"P = _ TT • (y _TT 

ZCF DCZ C LCZ 

F ZPF = _F DPZ* a P _F LPZ 
F ZSF = -F DSz’ a S -F LSZ 
F ZMF = _F DMF* a M~ F LMZ 



(3-19) 

(3-20) 

(3-21) 

(3-22) 



Equation (2-44) shows the general expression for hydro- 
dynamic forces, and from Ref. 3 the linear approximation of 
hydrodynamic coefficients are: 



CL cf =5.15 


a c 


+ 


0.036 


S F + 0.13 


CL pF =4.10 


a P 


+ 


0. 018 


6 P 


CL gF =4.10 


a s 


+ 


0. 018 


6 S 


CL irc , = 5.15 
MF 


a M 


+ 


0. 036 


6 a 
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Let 



C DF = °- 012 + 0.0015 < 5 p 

C DA = °’ 02 + °’ 0028 6 a 
C DP = 0,02 + °- 0028 6 p 
C DS = 0,02 + °* 0028 6 S 

0.5*p*U 2 *A CF = D i 

0.5*p*U 2 * A pp = 0 . 5 * P * U 2 *A gF = D 2 



After dropping the nonlinear terms 

F ZCF = -( 5 - 162D i )w'+(5. 162Di )Q 1 -D 1 (0. OOISA^O. 036)S F 
— D , ( 5 . 162 A j +0 . 13 ) 

F ZPF = -( 4 - 1 2D 2 )w'+(4.12D 2 )Q 2 -4.12D 2 A 2 
F ZSF = "(4.12D 2 )W'+(4.12D 2 )Q 2 -4.12D 2 A 2 

F ZMF = -( 5 - 17D 3 ) w ' + (5.17D3)Q 2 -D 3 (0.0028A 2 +0.036)6 A -5.17D 3 A 2 

w 

F z = -(5.162D 1 +8.24D 2 +5.17D 3 )§ +( 214 . 69DJ -139 . 92D 2 -87 . 79D 3 ) 

-D! ( 0 . 001 5Aj +0 . 036 )6 p -D 3 ( 0 . 0028A 2 +0 . 036 )6 a 

-[D! (5.162Aj+0.13)+8.24D 2 A 2 +5.17D 3 A 2 ] 

" = (5. 162D 1 +8. 24D,+5. 17D, ) * (214.6901-139.920,-87.790,) 

mU “ + mU 

1TT ,A. F G „ D t ( 0 . 001 5 A j +0 . 036 ) ,, D , ( 0 . 0028A, +0 . 036 ) , 

1 m m F m A 

[D, (5. 162A, +0. 13)+8 . 24D 0 A ? +5. 17D, A, ] 



m 



( 3 - 23 ) 



clUo 



M m fwd m aft' ! ' f xcf l zcf +f xsf l zsf +f xmf l zmf +f xcs l zcs +f xps l zps 

+f xss l zss +t x l zt 



M = F L 
FWD ZCF XCF 



AFT 



( F +F +F ~)h 
K ZPF ZSF ZMF ' XMF 



= [ (214 . 69D 3 )W -(214.69D 1 )Q 1 +(0.212D 1 A 1 +1.497D 1 )6 
+ ( 214-690^! +5. 41Di )] + [-( 139 . 92D Z )W +(139. 92D 2 )Q 2 
-139.92D 2 A 2 -(87.79D 3 )W ' +( 87 . 79D 3 )Q 2 - ( 0 . 048D 3 A 2 +0 . 611D 3 )6 
-87 . 79D 3 A 2 ] +[ (281. 1D 2 A 2 )w' - ( 281 . 1D 2 A 2 )Q 2 +(140. 55D 2 a! 
-0.686D 2 )]+[ (159. 65DiAi+2. 02Di )W ' - ( 159 . 65Di Aj +2 . 02Di )Qi 
+ (0. 558Di Ai-0. 023Di ) 6 p +(79 . 83Di A? +2 . 02AiDi -0 . 186 ) ] 

+ t (176.5D 3 A 2 )W-(176.5D 3 A 2 )Q 2 +(0.62D 3 A 2 -0.048D 3 )6 



Q 



+(88 . 27D 3 A 2 -0. 343)]-14.29U 2 +17. 7T 



X 



[D, (216. 71+159. 65A,)-D 2 (139. 92-281 . 1A, )-D ,( 87 . 79-179 . 5A, ) l w 

I U . " 

X YY 

(Pi (8928. 96+6639. 84A, )+D, (2375 . 84-4773 . 08A, )+D 3 (1490. 67 

J YY U 

-2996. 97A ? )] + D, ( 1 . 474+0 . 77A 3 ) , _ D, (0. 659-0. 572A, ) fi 



YY 



YY 



+ tDi 5.41{Ai (216. 71+79. 83A 1 ) > -P 2 { 0 . 686+A 2 (139.92 



-140 . 548A 2 )}-D 3 A 2 (87. 79-88. 27A 2 )-0. 529-14 . 29U 2 +17 . 7T x ] 



YY 



(3-24) 
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From (3-23) and (3-24) 



e + c 2 e - CjZ = c 3 <5p, - c 4 <s A + Uj 
-c 6 e - c 7 e + z + c 5 z = -c 8 6 F -c 9 6 A - u 2 



By taking the Laplace Transform of the two equations 
above and rearranging into a matrix form 



' 

S(S+C 2 ) 


/ 

CO 

H 

0 

1 




e 




C 3 


-c 4 ' 


V 




'Uj ’ 


-(C 6 S+C 7 ) 

V 


S(S+C 5 ) 




z 


= 


— C 8 


1 

O 

<x> 

v - - 


Cb 

> 


+ 


-U 2 



The characteristic equation 

A = S[S{S 2 +(C2+C 5 }-CiC 6 S-CiC7] 



wnere : 



C, = 



[ D i (216.71 139. 65A! )-P 2 (139. 92-281. 1A ? )-P 3 (87 . 79-176 . 5A 2 ) ] 



I YY u 



„ _ [D 2 (8928. 96+6639. 84Aj )+D 2 (2375. 84-4773. 08A 2 ) 

- i u 

iY 

+D, (1490. 67-2996. 97A 2 ) ] 



P _ [D^l. 474+0. 77A,)1 

J 

i YY 



„ _ [D 3 (0. 659-0. 572A 2 )] 

- j 

YY 



Cs = 



_ [ 5 . 162Di +8.24D 2 +5.17D 3 ] 



mU 



„ _ [ (214.69Di-139.92D2-87.79D3 ) +U] 

C 6 
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/ 



c - F ° 

C7 - — - 
m 



_ [Di (0. 0015 A i+O. 036) ] 
'-'8 ~ 



m 



After inserting the value of parameters, a speed of 36 
knots was chosen to be used in this simulation to insure 
that cavitation did not exist, and it was assumed that there 
was no external disturbance. The transfer function matrix 
becomes : 



f 


X 




r s 




t 


> 




/ V. 


S(S+10. 791) 


-0. 063S 




e 




0.083 


-0.0394 




S F 


- (76. 33S+32 . 2 ) 

\ 


S( S+7 . 61 ) 

/ 




z 

i 




-1.04 


-1.215 




. 6 a. 



and the characteristic equation was 
A = S[S 3 +18.4S 2 +77. 31S-2. 03] 



= S( S-0. 0261 ) ( S+6. 55 ) ( S+ll .876) 



and 



0.083 -0.063S 

JL 

j? _ 

6 A =0 -1.04 S(S+7. 61) 

A 



0 . 083( S+6 . 82 ) 

(S+0. 0261) (S+6. 55)(S+11. 876) 




S( S+10. 731) 0.083 

(76 . 33S+32 . 2 ) 1.04 

A 



1. 04(S-0.495)(S+5. 19) 
S(S-0. 0261) (S+6. 55)(S+11.876) 
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-0.0394 -0.063S 

-1.215 S(S+7. 61) 

A 



-0. 0394(S+9. 55) 

(S-0. 0261)(S+6. 55)(S+11. 876) 




S(S+10 . 791 ) 0.0394 



-(76. 33S+32. 2) 1.215 

A 



1.215(S+0.079)(S+13.186) 
S(S-0. 0261 ) (S+6 . 55 ) ( S+ll . 876 ) 



Wg = -Usin0+Vcos0sin<}>+Wcos0cos<j) 

V, the velocity in Y-direction is equal to zero, and 
angle 0 is assumed small. 

Sin 0=0, Cos 0=1. 



W £ = -U*0 + W 



- h cg ■ - Z E - -<!- I‘ e > 



The height at the sensor, 
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The four transfer functions of the craft are: 

0. 083(S+6. 82) 

(S-0.0261)(S+6. 55)(S+11. 876) 

6. 68(S+0. 72)(S+6. 52) 

S(S-0. 0261)(S+6. 55)(S+11.876) 

-0. 0394(S+9. 55) 

(S-0. 0261)(S+6. 55)(S+11.876) 

-1.464(S+2. 72)(S+5. 348) 

S ( S-0. 0261) (S+6. 55)(S+11.876) 

The characteristic equation shows that this system has 
one pole jin the right half plane, and from the Root Locus 
plots three of the plant transfer functions have roots in 
right half plane, only the effect of forward foil to the 
height sensor has complex roots in the left halp plane but 
very close to the imaginary axis, which was only marginally 
stable. From this point, the stability of the craft without 
automatic control system was questioned, and led to the con- 
clusion that the craft must have an automatic control system. 

B. ROOT LOCUS DESIGN OF FORWARD FOIL AUTOMATIC CONTROL 

Considering only the forward foil as a control surface, 
the aft foil was fixed at the angle of 0.074269 radian [Ref. 
1]. By closing the height loop and the pitch angle loop as 
in Figure 14 and increasing gain K of the system to 5, from 
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Figure 10 the roots of the height loop moved further into 
the left half plane and became repeated real roots at -1.6. 
This value of gain K pushes the root of the pitch angle loop 
(Figure 11) just across the imaginary axis to the left half 
plane. The system becomes stable. The response to the 
change in height command of 2.0 feet was shown in Figure 15. 
The craft went to steady state in 5.5 sec. 




FIGURE 14 

BLOCK DIAGRAM SHOWING FORWARD FOIL AUTOMATIC CONTROL 

C. ROOT LOCUS DESIGN OF AFT FOIL AUTOMATIC CONTROL 

In like manner, if the aft foil is used as a control 
surface and the forward foil is fixed at the angle of 
0.03711 radian, the system is unstable by itself, because 
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X-Scale 6.0 
Y-Scale 1.0 



FIGURE 15 

STEP RESPONSE OF THE CRAFT 
(FORWARD FOIL ONLY) 
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the open loop roots of both height and pitch angle were in 
the right half plane. From the Root Locus plots Figures 12, 
13, positive gain will have the roots further into the right 
half plane. After closing the loops and using gain K equal 
to -47, the height loop had complex roots at -3±j3. The 
right half plane root of the pitch angle loop moved toward 
the left half plane, across the imaginary axis to the point 
-0.2. The system becomes stable. (See Figure 16.) 




FIGURE 16 

BLOCK DIAGRAM SHOWING AFT FOIL AUTOMATIC CONTROL 

Figure 17 shows the response of the craft to the change 
in height command of two feet. Because the flap angle was 
limited at ±5°, and the gain of the aft foil was too high, 
in the first five seconds the flap exceeded its limit and 
the flap angle was set to 5° which had the same result as 
reducing the gain of the system. Then in the first five 
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Y-Scale 1.0 

FIGURE 17 

STEP RESPONSE OF THE CRAFT 
(AFT. FOIL ONLY, INPUT = 2.0) 
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seconds the craft behaved like a very lightly damped system, 
but after five seconds the .craft.; became more stable and 
reached steady state in eight seconds. Figure 18 shows the 
response of the craft when the change in height command was 
0.5 foot. In the first 1.0 sec. the flap angle still ex- 
ceeded the limit, but the craft went to steady state in 3.0 
seconds. 

Parts B and C have shown that the stability of the craft 
in pitch-heave mode can be achieved by using either forward 
or aft foil as a control surface with the other foil fixed. 
But to make the craft completely stable in six degrees of 
freedom both forward and aft foil have to be used as control 
surfaces. The dynamics of the system will produce a cross 
coupling .effect between the two inputs, the systems now have 
multiple inputs and multiple outputs, which is a classical 
multivariable system and will be studied in the next chapter. 
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FIGURE 18 

STEP RESPONSE OF THE CRAFT 
(AFT. FOIL ONLY, INPUT = 0.5) 
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V. MULTIVARIABLE SYSTEM 



A. BASIC CONCEPT 

A system which has multiple inputs-mult iple outputs is 
much more complicated to design than a single input-single 
output system, because a change in one input may have an ef- 
fect not only on the corresponding output but also on the 
other outputs. In designing this kind of system, special 
techniques are needed. Ref. 9 has all of the theories and 
procedures to work on this kind of problem, and this paper 
will follow the procedures in that reference. 

The idea was that, the plant has a nxm transfer function 
matrix and was compensated by a iPn diagonal matrix compen- 
sator with unity feedback. The configuration of the compen- 
sated system is shown in Figure 19. 

Consider the hydrofoil craft in pitch-heave made only, 
and not coupled to the other modes. The craft will have a 
2x2 transfer function matrix. 



g Pll g P 12 

G = 

P . 

P* g 

P21 P22J 

If the diagonal compensator matrix 



G 

c 



g 



c 1 1 
0 



0 



g 



C 22 



(4-1) 



(4-2) 
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is used, the system configuration becomes 




FIGURE 19 

BLOCK DIAGRAM SHOWING A MULTIVARIABLE SYSTEM 
WITH DIAGONAL MATRIX CASCADE COMPENSATOR 



G G 
P 



c 



g pi l g C] 1 
V g P2 l g Cx 1 



g p 12 g c 22 

g P22 g C 22 



(4-3) 



1 . The Characteristic Equation 

The stability of a multivariable system is determined 
by the zeros of Ni(s) and N(s), which are defined as follows: 



57(§y - det U+G p (s)G c ( S )] 



N(s ) A 



A c (s)-A p (s) 



D i (s ) 



(4-4) 

(4-5) 



and A c (s), A p (s) represent the characteristic polynomials of 

the rational transfer function matrices G (s) and G (s), 

c v ' p ’ 

respectively . 

If 

1) Cancellations are selected systematically, as in 
Ref. 9, 
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2) Poles of G c are carefully selected, 
then (4-4) alone determines the stability of the system. 
Thus 



det ( I+G G ) = 0 

pc 



will be referred to as the characteristic equation of multi- 
variable system. 

2. Connection Between Single- Loop and Multivariable 
Cases for 2x2 Plant G , and 2x2 Diagonal G c 



det ( I+G G) = 1+g g +g g +(det • G )g g 

P C Pi 1 C 1 1 P2 2 C 2 2 P C 1 1 C 2 2 



g +(det • G )g 

■ 1+ it ; ,; 



Pi i c i i 



(4-6) 



- d+g D S. )(1+G en g. ) 

Pi 1 C 1 1 e Q C 2 2 



(4-7) 



where 



g n + ( det ‘ G n )g r 
G = , -P. 2 z 2 g±-L 

eq 1+g g 

Pi l c i l 



det • G 

g 1+ E 

P2 2 g„ 

P2 2 



1 1 



1+g n S c 
Pi l c l l 



The roots of the characteristic equation are the 

zeros of the rational function (1+g g )(1+G g ). 

Pi i Ci i eq C 22 

From (4-6), can be seen that the roots of 



det *G 

1 + 2 g 

g n Cl1 

P2 2 



0 



and 



(4-8) 



1 + 



g Pi l g Ci 1 



= 0 



(4-9) 
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constitute part of the zeros and poles of G , respectively 
det *G 

and £ and g are known functions. 

g P *2 P1 ‘ 

Therefore, any single-loop method can be used in de- 
signing g to place the roots of (4-8) and (4-9) at desir- 
c i i 

able locations, and also has toform a reasonably good 

pole-zero pattern for G . After g is designed, all 

eq ci i 

poles and zeros of G are known and g can be designed 
* eq & c 22 

from 



1+G g =0 

eq & c 2 2 



(4-10) 



B. ROOT LOCUS DESIGN FOR PITCH-HEAVE MODE 
From the craft transfer function matrix 



G = 
P 



-h. 



g 



= ^S = 6.68(S+0. 72)(5+6.52) 



Pi l 



S(S-0. 0261)(S+6. 55)(S+11. 876) 



g 



= = -1.464(S+2.72)(S+5. 348) 



Pi 2 



S(S-0. 0261)(S+6.55)(S+11. 876) 



g 



0. 083(S+6 . 82) 



P 2 1 



S p (S-0.0261)(S+6.55)(S+11.876) 



g 



•0. 0394(S+9. 55) 



p 2 2 <5 a (S-0. 0261)(S+6. 55)(S+11.87G) 
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det • G = g g -g g_ 

P P 11 P 22 P 12 P 21 



-0.141(S-0. 02325)(S+6. 51)(S+12. 02) 
S(S-0.0261) 2 ( S+6 . 55) 2 (S+11.876) 2 



det*G 



^ _ 3.579(S-0.02325)(S+6.51)(S+12.02) 



g 



P 22 



S(S-0.0261)(S+6.55)(S+11.876) 



The design procedure starting from g , by plotting the 

Ci 1 

root locus of equations (4-8) and (4-9) as in Figures 20 and 
21, from Figure 20 can be seen that, as the gain was in- 
creased the root in the right half plane moved toward the 
left half plane, met pole at the origin, broke away and 
circled around the zero at -0.72, with the value of the gain 
equal to 5. Two dominant roots were at -1.25 and -2.2, and 
other two roots at -6.75 and -11.45 and the effect of this 
gain pushes the right half plane root of equation (4-8) to 
the left and was cancelled by the zero at +0.02325 as in 
Figure 21, the pole at the origin moved further into the 
left half plane to the point -2.3 and the other roots were 
in the left half plane at -6.8 ± j0.8 and -11.99. Therefore, 

by using the gain equal to 5 as g , the height loop was 

c 1 1 

stabilized and, this also gave a good pole-zero pattern for 
G 



eq 



Equations (6-30) - (6-35) of Ref. 9 show the procedures 



to find poles and zeros of G , which were listed below 

eq 



Gi = 1+g g 

Pi 1 Ci 1 



G 2 = g^ +(det • G )g 

P2 2 P C 1 1 



(4-11) 



(4-12) 
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r - G z 
eq Gx 



P = Z r +P p -(P r 0 P p ) 
U e q Ui G 2 G 2 Gi 



(4-13) 

(4-14) 



Z G ea = Z G 2 +P G,- CP G 2 P Gl ) 
«q 



p = P 

Gl g Pii g cn 



= P +P 

g Pn g cn 



(4-15) 

(4-16) 



Z c = Zi+g g 

^1 P 1 i c i i 



(4-17) 



j g 2 



det • G 

= Zl + £ 



g. 



g n "Cll 

P 2 2 



(4-18) 



And from the system transfer function matrix 

= P- P = P +P = {0,0. 0261,-6. 55,-11. 876}+P 

Gl g n g c g n g r g e 

Pii c ii Pii c ii °i l 



Z r = Z x +g g 
Gi Pi l Cii 



P = P +P(det-G )g -[P 0 P(det*G )g ] 

U 2 g» P Ci i g P C x i 

P 2 2 P 2 2 

= { 0.0261, -6. 55, -11. 876} + {0, (0.0261) 2 , (-6. 55) 2 , (-11. 876 ) 2 } 
+P - [ { 0 . 0261 ,-6. 55,-11. 876}D { 0 , ( 0 . 0261 ) 2 , (-6 . 55 ) 2 , 

g Ci i 

(-H-876) 2 }+g ] 

c i i 

= {0, (0.0261) 3 , (-6.55) 3 , (-11.876) 3 +P -[{(0.0261), 

g cu 



(-6.55), (-11.876)}] 
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= {0, (0.0261) 2 , (-6. 55) 2 , (-11 . 876 ) 2 }+P 



g 



Cl l 



det • G 

z g 2 ■ z ' + -g— 11 « 



P2 2 



Cl 1 



P r 0 P r = (0,0.0261,-6. 55,-11. 876}+P = P r 

Gi G 2 g G 

Cl l 1 



P = Z,+g g +{0, (0 . 0261) 2 ,(-6.55) 2 ,(-11.S76) z }+P 
°eq p “ Cl > s o,, 



-{0,0.0261,-6. 55,-11. 3761-P 



g 



Cl l 



= Zi+g„ g„ +{(0.0261), (-6.55), (-11.876)} 

P 1 1 Cl 1 



det • G 

Zp = Z r — Z i + — — g 

G G 2 g Ci 1 

eq z & p 2 2 



G has three known poles at 0.0261, -6.55, -11.876; the 
other poles were the roots of equation (4-9) 



l+g n g_ =0 
Pi 1 c i 1 



The zeros of G were the roots of equation (4-8) 



det • G 

1+ £ g =0 

g P22 Cl1 



and therefore 



K . (S+2.3)(S+6.8±j0.8)(S+11.99) 
_C£L 



eq (S-0. 0261)(S+1.25)(S+2. 2 ) (S+6 . 55) (S+6 . 75 ) (S+ll . 45) (S+ll . 876 ) 



and from 
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g +(det*G )g 
G = P 22 P cn 

eq l+g„_ . g. 



pi 1 Cl i 



0. 0394 ( S+9 . 55 ) 0. 144 (S-0. 02325 ) ( S+6 . 51 ) ( S+12 . 02 ) 

(S-0. 0 26 1)( S+6. 55)( S+11.876) - S (S-0. 0261 ) 2 (S+6. 55 ) 2 (S+11.876) 2 

1 6. 68(S+0. 72) (S+6. 52) 

S(S-0. 0261) (S+6. 55) (S+11.876) 



-0. 0394 ( S+Z i )(S+Z 2 ) 
(S+Pi )(S+P 2 ) 



K = -0.0394 
eq 



(4-19) 



a -0 . 0394 ( S+2 . 3 ) ( S+6 . 8±.j 0. 8 ) ( S+ll . 99 ) 

eq ( S-0. 0261 )(S+1. 25) (S+2. 2) (S+6. 55) (S+6. 75) (S+ll. 65) (S+ll. 876) 



Equations (4-10) and (4-19) were the equations used in 
designing the compensator g 

P 2 2 

Figure 22 is the root locus plot for designing g 

P 2 2 

The positive gain will push the pole at 0.0261 further into 

the right half plane; therefore positive gain cannot be used 

to compensate this system. Negative gain which can make the 

gain of the system equal to 3.881 will set dominant roots at 

-0.7±j0.2, which is a pretty good position. Therefore, 

using the gain of the system equal to 3.881, we can make 

system stable and the pitch angle loop has good damping. 

The compensator g thus becomes a pure gain of value 

P 2 2 

K2 K =3.881 
eq 

3.881 



K2 



-0.0396 
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The configuration of the compensated system is shown in 



Figure 23. 




FIGURE 23 

BLOCK DIAGRAM SHOWING THE COMPENSATED SYSTEM 

Referring to Figure 20, the dominant roots were at -1.25 
and -2.2 which gave = 1.658, and = 1.04. Figure 24 
shows the response of the craft when the height above the 
water has changed from 9.5 feet to 11.5 feet, and Figure 25, 
when it changes from 11.5 feet to 9.5 feet. The response of 
the craft was the same in both cases and followed very 
closely to the response of the best second order approxima- 
tion . 

Figure 26 is the response of the pitch angle from 2° bow 
up back to 0° reference and Figure 27 is the response from 
2° stern up back to 0° reference. The response in both 
cases was the same but faster than the response of the best 
second order approximation. 
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FIGURE 24 

STEP RESPONSE OF HEIGHT (2 FT. UP) 
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Y-Scale 1.0 



FIGURE 25 

STEP RESPONSE OF HEIGHT (2 FT. DOWN) 
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Y-Scale 1.0 

FIGURE 26 

STEP RESPONSE OF PITCH ANGLE (BOW UP) 
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Y-Scale 1.0 

FIGURE 27 

STEP RESPONSE OF PITCH ANGLE (STERN UP) 
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C. ROOT LOCUS DESIGN FOR ROLL MODE 



After neglecting the small terms, equation (2-4) becomes 




And from equation (2-41) 

L = f zpf l ypf +f zsf l ysf -f yps l zps _f yss l zss -f ycs l zcs 

Assume that the roll mode was decoupled from the other 
modes 

L = f zpf l ypf +f zsf l ysf 

= [-D 2 (0. 0028A 2 +0. 018 )6 p ] L ypF + [ -D 2 ( 0 . 0028A 2 +0 . 018 ) 6 g ] L ygF 
S p = -« s 

L = -29959. 46 g +(-29959. 46 g ) 

= -59918. 86 g 

P = i = -0. 16g 

Taking the Laplace Transform, the transfer function is 
simply 




Figure 28 is the root locus plot after compensation 
using a lead network which has zero at -1.5 and pole at -15. 
With the gain of -82, the roots were at -3.5±j0.75, which 
gave the value of £ = 0.96 and = 3.6. Figure 29 is the 
step response roll from 2° port up back to 0° reference and 
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Figure 30 from 2° starboard up to 0° reference. The same 
response is obtained in both cases. The block diagram of 
the automatic control of this system is shown in Figure 31. 
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FIGURE 29 

STEP RESPONSE OF ROLL ANGLE (PORT UP) 
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FIGURE 30 

STEP RESPONSE OF ROLL ANGLE (STARBOARD UP) 
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BLOCK DIAGRAM OF CONTROL SYSTEM 



VI. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY 



The transfer functions and root locus study has shown 
that this system is unstable by itself because of the open 
loop pole in the right half plane. From studying the effect 
of each foil on the system, it has been shown that stabiliza- 
tion can be achieved by simply closing each loop and using 
the appropriate gain. However, that will stabilize only the 
effect of one foil when the other is fixed. 

To stabilize the systems when both forward and aft foils 
are used, the total system has been treated as a multivari- 
able system. To reduce the work in designing the automatic 
control system, some assumptions had to be. made, such as no 
coupling between the pitch-heave modes, roll mode and yaw 
mode. Thus, the system will have only a 2x2 transfer func- 
tion matrix for the pitch-heave mode, and single loop for 
the roll mode. The step response study has shown that the 
responses of the compensated system followed closely the re- 
sponse of the best second order approximation. The system 
was made stable in both pitch-heave and roll mode. However 
this result cannot insure that the system is completely 
stable in six degrees of freedom. 

Recommendations for further study are: 

1) Study of yaw motion and the craft turning character- 
istic 

2) Wind effects 
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3) Produce a simulation for the transition period be- 
tween hullborne and foilborne 

4) Take the dynamic characteristic of the sensors into 
account 

5) Study the response of the craft in a sinusoidal sea 
and in a random sea. 



I 
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APPENDIX A 



SUMMARY OF COMPLETE EQUATIONS 
EQUATIONS FOR ACCELERATION IN BODY AXIS 



U = ^F+RV-QW 

m A 


( A-l) 


V = — F+PW-RU 
m Y 


(A-2) 


W = ^F+QU-PV 
m ia 


(A-3) 


P = I i-[L-QR(I zz -I yy )+(R+QP)I xz ] 

XX 


(A-4) 


R = Y^[N-PQ(I yy -I xx )-(QR-P)I xz ] 


(A- 5) 


VELOCITY IN BODY AXES 
U = / Udt 


(A-6) 


V = / Vdt 


(A~7) 


W = / Wdt 


(A-8) 


P = / Pdt 


( A-9 ) 


Q = / Qdt 


( A-10) 


R = / Rdt 

TRANSFORMATION OF VELOCITIES FROM BODY TO EARTH AXES 


( A-ll ) 



U v =Ucos0cosi{j+V(cost{isin0 sincfi-simjjcostf) ) +W( cos^sin0cos<t> + sinijjsin<j) ) 
h (A-12) 

V v =Ucos0simfi+V(cos^cos<j>+sirnJisin0sin<j) )+W(sin<jJsin0cos4>-cosijjsin(b) 
* (A-13) 

Wg=-Usin0+Vcos0sin(})+Wcos0cos(J) ( A- 14 ) 
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1 

} = P+4»sin9 ( A-15 ) 

9 = Qcos<f)-Rsin(() (A-16) 

ip = (Qsin<j>+Rcos<j> )cos9+(|-P)sin9 (A-17) 

POSITION IN EARTH AXES 

X E = / U E dt ( A-18 ) 

Y E = / V E dt ( A-19 ) 

Z E = / V/ £ dt ( A-20) 

<p = / 4>dt ( A-21 ) 

9 = / 9dt ( A-22 ) 

ip = / if>dt (A-23) 

EXPANSION OF FORCE AND MOMENT EQUATIONS 

F X = E F XIF + E F XIS + mg X +T X ( A-24 ) 

F y = E F YIS + mg Y ( A-25 ) 

F Z = E F ZIF + mg Z ( A-26 ) 

L = ( f 2 Ij Y^PF + ^ F Z Ij Y^SF _ ^ F Y L Z^PS“^ F Y L Z^SS _ ^ F Y Ij Z' ) CS ( A ~ 27 ) 

M = 2 <Wif + i: ( f xVif +2 <Wis +t x l zt (A - 28) 

N = E (F y L x ) IS ( A-29) 

EXPANSION OF GRAVITY TERMS 

mg^ = -mg sin9 (A-30) 

mg Y = mg cos9sin<{> (A-31) 

mg^ = mg cos9sin(f> (A-32) 
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FOIL VELOCITY COMPONENTS IN BODY AXES 



Center Foil: 

u c = U+L ZCF* q 

v c = v - l zcf- p+l xcf- q 

w c = W - L XCP ,Q 



Port Foil: 



U T 



V T 



W T 



u+l zpf*^ ■Spf* 11 
v_l zpf’ p+l xpf‘ r 

W+L YPF ,P - L XPF* Q 



Starboard Foil: 

u s ' u+l zsf'«- l ysf- s 
v s = v " l zsf' p+l xsf' r 

W S - W+L YSF- P - L XSF-« 



Mid Foil: 



U M = U+L 



'XMF 



Q 



= V-L • P+L 
M ZMF XMF 

= W-L 

M XMF 



( A-33) 
( A-34 ) 
(A-35) 



( A-36 ) 
( A-37 ) 
( A-38 ) 



( A-39) 
( A-40) 
( A-41) 



( A-42 ) 
( A-43 ) 
( A-44 ) 
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TRANSFORMATION OF WATER PARTICLE ORBITAL VELOCITY FROM EARTH 



TO BODY AXES 

U WI = u EW i cos ^ cose+V EWI Sin ^ COS0-W EWI Sin9 ( A-45 ) 

V WI = U EWI (cosi|;sin0sin<{)-sim(jcos({) )+V WEI (sin^sin0sinc})+cos^cos4) ) 

+ ^EWI COs9s '*‘ n< ^ (A-46) 

^WI = U E ^j(cosiJ;sin0cos(|)+sim(jsin({) )+V^ E j (sim|)sin0coscf)-cosi|'sin(f) ) 

+ ^EWI COs9cOS ^ (A-47) 

Each of these equations must be repeated for each foil. 



RELATIVE VELOCITY COMPONENTS 



U RI U I U WI 



V -= v - V 
V RI I WI 



w - \v - W 
RI I WI 



One set of equations for each foil 



( A-48 ) 



( A-4S ) 



(A-50) 



ANGLES OF ATTACK AND SIDE SLIP 



a I arctan W RI /U RI + a x_ FIX ED 



8 t = arcsin 



RI 



/u 2 +v 2 +w 2 
v RI RI RI 



(A-51)* 

(A-52)* 



* One set of equations for each foil. Angles of attack and 
side slip are calculated in radians. 
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TOTAL RELATIVE VELOCITY AT A PARTICULAR FOIL OR STRUT 



V = /U* + V* + W* ( A-53 ) 

HYDRODYNAMIC FORCES IN WATER AXES 

(Each foil equation adn strut equation must be repeated for 
foil and strut, respectively) 

Foil Lift: 

F l = -|-pV 2 A F C F (A- 54) 

Foil Drag: 

r DF = iP y2A F C F (A_55) 

Strut Drag: 

P DS ^ V2A S C DS (A " 56) 

I 

Strut Side Force: 

Fg = ipV 2 AgCg ( A-57 ) 

TRANSFORMATION OF FOIL FORCES FROM WATER AXES TO BODY AXES 



XIF 


“ F IDF COSaCOS ^ +F LIF sina 


( A-58 ) 


YIF 


- F DIF sinB 


( A-59) 


ZIF 


-F DIF SinaCOS ^ -F LIF COSa 


( A-60) 
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APPENDIX B 



PLOTS OF HYDRODYNAMIC COEFFICIENTS 



Fwd. Foil Lift Coefficient vs. Angle of Attack 




Y-Scale 0.4 Unit/inch 
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Mid Foil Lift Coefficient vs. Angle of Attack 




Y-Scalo 0*4 Unit/inch 
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Aft. Foil Lift Coefficient vs. Angle of Attack 




Y-Scale 0 o 4 Unit/inch 



91 



a 00 Xt-OO 6-00 300 10-00 ilOO 

(f 10 ) 



Forward Foil Drag Coefficient vs. Angle of Attack 
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// EXEC OSL » P ARM.C = 'NOS OURC E , NOMA P 1 
//DSL. INPUT DO * 

D COMMON MATRI X ( 53 9 ) 

CONTRL FINTIM=10.0,DELT=0.025,DELS=0.05 
INTEG RKSFX 

INTGER IFLAG , I FLAG I , NPLOT »NUM»NCUR 
CONST NPLOT=3 ,NCUft=l , IFLAG=0 

CONST WEIGHT = 26252 7, IXX=60I200, IYY=4286000, . . . 

IZZ=41 74000, G=32. 2 

CONST ACFP=32.7,ACFS=32.7,APF=36.8,ASF=36.8,AMF=76.0, ... 
PI=3. 141593 

CONST LXCF=41 . 59 , LXCS =4 1 . 59 , LX PF = -1 6 . 98 , LXPS=- 16 . 93 , ... 
LXSF=-lo.98,LXSS=-16.98,LXMF=-16. 98, LYCP=-5.0, LYCS=5.0, . . . 
LYPF=-12. 42 ,LYPS=-8.0,LYSF=12.42,LYSS=8.0,LZCF=15 .5, . . . 
LZCS=13.25,LZPF=17. 14, LZ P S= 1 4. 0 , L Z SF= 17. 14 , LZ S S = 1 4. 0 , . . . 

LZMF = 17 .14,LXT=-16«9,LZT=17 . 7 , LX H=68 . 0 
CONST K1=5.0,K2=— 98.0, K3 = -82.0 
STORAG IC ( 2 ) , NU ( 1) ,0EN(3) 

TABLE I C ( 1-2) =2*0.0,NU( 1 ) = 3. 99, DENI 1-3) =1.0, 4. 0,3. 99 
PARAM U ZERO= 36 . 0 

PREPAR .1,FX,FY,FZ,L,M,N,U,V,W,ZE,P,Q,R, THETA , UDOT , WDOT , . . . 
QDOT, ALFCP, ALFCS,ZE1 , PHI , PH 1 1 , THE TA3 , X , Y 
GRAPH TIME,THETA3,X,Y 
PRPLOT 

LABEL STEP RESPONSE OF PITCH ANGLE 

PRINT 0.5, TIME, ZE, PHI ,P , PDOT , TH E T A , EL VP , EL VS , F LAP , ARS , ARP 

nC 

INITIAL 

* 

DENSE = 2 .0 

UINIT=UZER0*2000. 0*3. 0*1.0/3600.0 

U=UINIT 

T X=3960 . 0 

A LM=0. 074269 

ALS=0. 074269 

ALP=0. 074269 

ALF=0. 037111 

PREF=0 • 0 

TREF=0 . 0 

THETA2=-2 . 0*P 1/180.0 

PHI 2=0 • 0 

MASS= WE IGHT/G 

CONV=P I /180.0 

FL AP=0 . 0 

E LV=0 « 0 

ARP=0 . 0 

ARS=0 . 0 

HREF= 1 1 .5 

HREF2=- 11.5 

FLAG=0. 0 

I FLAG 1 = 0 

EPS=25. 0 

jjc 

DERIVATIVE 

* 

NOSORT 

* ANGLES OF ATTACK AND SIDE SLIP 

* 

VRC=V-LZCF*P+LXCF*R 

VRS=V-lZSF*P+LXSF*R 

VRM=V 

VRP=V-LZPF*P+LXPF*R 

WRP=W+LYPF*P~LXPF*Q 

WRS=W+LYSF*P-LXSF*Q 

WRM=W— LXMF *Q 

WRC=W-LXCF*Q 

ALFC=WRC/U+ALF 

ALFAP=WRP/U+ALP 

ALFAS=WRS/U+ALS 

ALFAM=WRM/U+ALM 

8ETAC=VRC/U 
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BETA P=V RP/U 
BETAS=VRS/U 

❖ 

* CALCULATION OF EULER ANGLES 

* 

PHI DI=P 

THETD1=Q 

PSID1=R 

* 

SORT 

THETA=INTGRL(THETA2»THETDI) 

PHI=INTGRL(PHI2, PH 1 01 ) 

PSI=INTGRL(0.0,PSIDI) 

WE=-U*S INI THETA) +V*COS ( THET A ) * S I N ( PH I ) + W *COS< THETA) . .. 
*COS( PHI) 

ZE=INTGRL(HREF2»WE) 

ZE1=ZE-HREF2 

£ 

NOSORT 

* 

* SUBMERGENCE TERMS 

# 

SECP=ZE-LXCF*THETA+LYCP*SIN ( PHI ) +LZCF*C OS ( PH I ) 
SECS=ZE-LXCF*THETA + LYCS-SI!M ( PHI ) +LZCF*CGS ( PHI ) 
SEPF=ZE-LXPF*THETA+LYPF*SIN (PHI ) +LZPF#COS( PHI ) 
SESF=ZE-LXSF*THETA+LYSF*SINIPHI ) +LZSF*COS ( PHI ) 
SEMF=ZE-LXMF*THETA+LZMF^COS(PH I ) 

$ 

* HYDRODYNAMIC COEFFICIENTS 

* 

ELV1=ELV*C0NV 
E LVD=LI MIT(-5.0,5.0,ELV1) 

FLAD=FLAP*CONV 

FDEG=LI MIT (-5.0»5.0, FLAD) 

ARPD=ARP*CONV 

APP1 = L J MIT (-5.0, 5.0, ARPD) 

ARSD=ARS*CO:MV 

ARS1=LIMIT(-5.0,5.0,ARSD) 

CALL INTERPI ALFC, ELV D, COEFC , 2 ) 

CALL INT£RPIALFAP,ARP1 ,CGEFPX,3 ) 

CALL INTERPI ALF AS , ARS1 , COEF SX , 3 ) 

CALL INTERPI ALFAM , FLAP , COEFMX , 1 ) 

CALL INTRP1(BETAC,S£CF,CFDS,2, THETA, THETA1 ) 

CALL INTRPK BET AP t SEP F,CPDS, I» THETA, THET A1 ) 

CALL INTRPK BETAS, SESFfCSDS, 1 , THETA, THET A1 ) 

* FORCES IN THE WATER AXES 

# 

FGRAX=-MASS*G*SIN( THETA) 

FGRAY=MASS*G*SIN(PHI) 

FGRAZ=MASS*G*COS( PHI ) 

CCS=- 2. 15*BETAC 
ACCS=CCS*4.25*SECF 
F SC = 0. 5*DENSE*IU**2)*ACCS 
FDSC=0. 5*DENSE*I U**2) *CFDS 
PCS=-2. 15*BET AP 
AFC$= PC S#4 .25*SE PF 
FPS=0.5*DENSE*IU**2)*APCS 
FDPS=0.5*DENSE*(U**2)*CPDS 
SCS=- 2. 1 5* BET AS 
ASCS=SCS*4.25*SESF 
F SS=0. 5*CENSE*(U**2)*ASCS 
FDSS=0 o 5*DENSE*( U**2) *CSDS 
C LCF = 5 . I5*ALFC + 0 .036*ELV+0. 13 
FLC=0. 5*DENSE*IU**2)*ACF*CLCF 
FDC=0. 5*DENSE *( U**2 ) *ACF*COEFC 
C LPF=4. 10*ALFAP+0.018*ARP 
FL P=0. 5*DENSE*I U**2) *APF*CLPF 
FDP=0. 5 ODENSE *( U**2 ) *APF*COEFPX 
CLSF=4. 10*ALFAS+0.018*ARS 
FLS=0.5*DENSE*(U**2)*ASF*CLSF 
F D S = 0 . 5 * D E N S E * ( U * * 2 ) * A S F * C 0 E F S X 
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# # # •K’C n •!{•*■«■ ■&**■«• 



CLMF=5. 1 5#ALF AM+O .036#FLAP 
FLM=0.5*DENSE*(U**2)*AMF*CLMF 
FDM=0. 5*DENSE *( U**2 ) *AMF*COEFMX 
FLCZ=0. 5*DENSE*(U**2)*ACF*CLCF 
FDCZ=0. 5*DENSE* ( U**2 ) *ACF*COEFC 
FLPZ=0. 5*DENSE*(U**2)*APF*CLPF 
FDPZ=0. 5*DENSE*( U**2) *APF*COEFPX 
FLSZ=0. 5* DENSE*! U**2 ) *ASF*C LSF 
FDSZ=0. 5*DENSE*< U**2) *ASF*CGEFSX 
FLMZ=0. 5*DENSE*< U**2) *AMF*CLMF 
FDMZ=0. 5*DENSE* ( U**2 ) *AMF*COEFMX 

TRANSFORMATION TO BODY AXES 

FXCF=-FDC+ALFC*FLC 

FXPF=-FDP+ALFAP*FLP 

FXSF=-FDS*ALFAS*FLS 

FXMF=-FDM+ALFAM*FLM 

FXCS=-FDSC-BETAC*FSC 

FXPS=-FDPS-BETAP*FPS 

FXSS=-FDS S-BE TA S*F SS 

FYCS=-BETAC*FDSC+FSC 

FYPS=-BETAP*FDPS+FPS 

FYSS=-B£TAS*FDSS*FSS 

F ZCF=- r FDC*ALFC _ FLCZ 

F ZPF=-F D PZ*ALFA P- F L P Z 

FZSF=-FDSZ*ALFAS-FLSZ 

F ZMF=-F DMZ*ALF AM-F LMZ 

LZCS=LZCF-0.4*SECF 

LZPS=LZPF~0.4*SEPF 

LZSS=LZSF-0.4*SESF 

FXONE=FXCF+FXPF+FXSF+FXMF+FXCS+FXPS+FXSS 

SUMMATION OF FORCES AND MOMENTS 

F‘X=FXONE+F GRAX+TX 
UDOT= ( 1 ./MASS)*FX-Q*W 
I F ( IFLAG1.EQ.O) UD0T=0„0 
FY=FYCS4-FYPS+FYSS + FGRAY 
FZ=FZCF+FZPF<-FZSF + FZMF + FGRAZ 

L=FZPF*LYPF+FZSF*LYSF-FYPS*LZPS-FYSS*LZSS~FYCS*LZCS 

MFWD=FZCF*LXCF 

M=-MFWD-FZPF*LXPF-FZSF*LXSF-FZMF*LXMF+. . . 

FXCF*LZ CF + FXPF*LZ PF + FX SF*LZ SF + FXMF*LZMF + FXC S. . . 

*LZCS+FXPS*LZPS*FXSS*LZSS+TX*LZT 

N=FYCS*LXCS+FYPS*LXPS+FYSS*LXSS 

VELOCITY INTEGRALS 



VDOT={ 1 . / M A S S ) * F Y - R*U 
WDOT= { 1 ./MASS) *FZ+Q*U 
PDOT= { 1 . 0/ I XX ) *L 
QDOT= ( 1 . 0 / I Y Y ) * M 
RDOT= ( 1 .0/IZZ)*N 
U=INTGRL(UINIT,UDOT) 

V=INTGRL( 0.0,VD0T ) 

W=INTGRL(0.0,KD0T) 

P=INTGRL(0.0,PD0T ) 

Q=INTGRL(0.0,QD0T) 

R=INTGRL(0.0,RD0T) 

AUTOMATIC CONTROL BLOCK 

❖ 

HACT=-ZE+LXH«THETA 

DELH=HREF-HACT 

DETH=TREF-THETA 

DEPH=PREFtPHI 

ELV=K1*DELFI 

FLAP=K2*DETH 
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ARS1 = L£DL AG( 0 .0 , 0 . 66 7 , 0 . 066 7 ,D E PH ) 

ARS=K3*ARS1 

ARP=-AR S 

THETA3=( THETA-THETA2) *180. 0/PI 
PHI 1=( PH I -PH 1 2) *1 80. 0/P I 
X1=STEP (0.0) 

X=2.0*X1 

Y=TRNFR (0 » 2» IC,NU,0EN,X) 

$ 

DYNAMIC 

$ 

SAMPLE 

CALL DRWGI 1 , 1 , T I ME , THETA3 ) 

CALL DRWG( 1 ,2, TIME, X) 

CALL DRWG< 1,3, TIME, Y) 

TERMINAL 

CALL ENDRW ( NPLOT ) 

END 

STOP 

FORTRAN 



SUBROUTINE SETUP 

THIS SUBROUTINE FORMS THE DATA POINTS COMPRISING 
THE CURVES OF FOIL AND STRUT ORAG COEFFICIENTS 
AND SIDE SLIP COEFFICIENTS INTO MATRACES. 



SUBROUTINE SETUP ( FLAG ) 

COMMON COEFMX, COE FCX, COEFOX ,CFDS1, CADS 1 
DIMENSION COEFCX (3,33) , COEFOX ( 3 , 3 3 ) , COEFMX ( 3 , 33 ) , 
1CFDS1 (lltll) tCAOSlf 11 til ) * INFOf 8) 

REAL INFO 

IC=0 

I FLAG=0 

I FLAG1=0 

KATRI X= 1 

J=0 

K=1 

IC0UNT=0 

2 READ(3,100,END=5Q) INFO 
100 FORMAT ( 8F10.5) 

I 1 = 1 

I F( IFLAG1 . EQ. 1 ) GO TO 51 
GO TO <4,5»6,7,8) » MATRIX 

4 CONTINUE 
DO 1 1=1,8 
J = J + 1 

I F( J.GT .3 ) K=K+1 
I F( J.GT .3 ) J = 1 
COEFOXt J,K)=INFO( I ) 

I COUNT = ICOUNT + 1 
IC=IC+1 

I F ( ICOUNT . EQ . 99 ) GO TO 3 
1 CONTINUE 
GO TO 2 

3 11 = 1 

I F ( 1 1 *L E. 8 ) 11=11 + 1 

J=0 

K= 1 

I COUNT= 0 
MATRI X= 2 

5 IF(11.GT.8) GO TO 2 
DO 9 1 = II , 8 

J = J + 1 

IFtJ.GT.3) K=K+ 1 
IF(J.GT.3) J=1 
COEFCX (J»K)=INFO( I ) 

I COUNT = ICOUNT + 1 
IC=IC+1 

I F ( ICOUNT. EQ. 99) GO TO 10 
9 CONTINUE 
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GO TO 2 
10 11=1 

IFU1.LE.8) 11 = 11 + 1 

J = 0 

K=1 

I COUNT=0 
MATRIX=3 

6 I F ( I 1 . GT. 8 ) GO TO 2 
DO 11 1=11,8 

J=J + 1 

IF(J.GT .3) K=K+ 1 
IF(J.GT.3) J=1 
COEFMX ( J,K)=INF0( I ) 

I COUNT= ICOUNT+1 
IC=IC+1 

I F C IC0UNT.Ea.99) GO TO 12 

11 CONTINUE 
GO TO 2 

12 11=1 

lFdl.LE.8) 11 = 11 + 1 

J=0 

K= 1 

I COUNT = 0 
MATRI X= 4 

7 I F ( II .GT.8 ) GO TO 2 
DO 13 1=11,8 

J = J + 1 

I F ( J .GT .11) K=K + 1 
IF(J.GT.ll) J=1 
CADS1 ( J » K ) = I NFO ( I ) 

I COUNT = ICOUNT + 1 
1 C=IC+1 

I F ( I COUNT .EQ.121) GO TO 15 

13 CONTINUE 
GO TO 2 

15 11=1 

lF(Il*LEo8) 11=11+1 

J = 0 

K=1 

I COUNT = 0 
MATR I X= 5 

8 I F ( II oGT.8) GO TO 2 
DO 14 1=11,8 

J = J + 1 

IF(JoGT.ll) K=K + 1 
IF(J.GT.ll) J=1 
CFDS1 ( J ,K) =1 NFO( II 
I COUNT = ICOUNT + 1 
IC=IC+1 

IF( ICOUNT. EQ.121) GO TO 16 

14 CONTINUE 
GO TO 2 

16 I FLAG1= 1 
GO TO 2 

50 I FLAG= 1 

IF( (IFLAG.EQ.l) .AND. I IFLAG1 .EQ.l) ) GO TO 52 

51 WRITE (6,104) ICOUNT , MATRI X , I C * 

104 FORMAT ( ' ' , ' JOB FA I LURE « , 1 X , 1 C OUNT= • , I 6 , IX , • MATR I X 
1I6,1X, • NG.PTS=' , 16) 

F LAG= 1 
RETURN 

52 WRITE(6,105) IC 

105 FORMAT ( * • , • JOB COMPL E TE 1 , 1 X, » PO I NTS ST0RED = , ,I6) 
RETURN 

END 

FORTRAN 



SUBROUTINE INTERP 



SUBROUTINE I NTER P ( AL FAC , F LA P , COEF , M ) 
COMMON COEFMX, COEFCX, COEFOX 
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DIMENSION COEFMXI 3,33) , COEF CXI 3 ,3 3) , COEF OXI3, 
DATA CONDEG/57. 295773/ 

ALFA=ALFAC*CONDEG 

I FI ALFA. GT. 10. 0) WRITE(6,100)ALFAC,M 
IFIALFA.GT .10.0) ALFA=10.0 
100 FORMAT! 'O' , ' ERROR. ALFA EXCEEDED L I M I TS. A LF A= 1 
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1 2X , ' M = ' ,111 
IFIALFA.LT. -6.0) 
I FIALFA.LT .-6.0) 
ALPHA=2 .0*ALFA 
IFIAL PHA.LT.O.O) 
IF(ALPHA.GT.O.O) 
I A=13 
WTA=0 • 0 
GO TO 3 
I Al=- AL PHA 
I A= 13- I A 1 
ALPHA1=-IA1 
WTA=ALPHA-ALPHA1 
GO TO 3 
I A1=ALPHA 
IA=13+IA1 
ALPHA 1= I A1 
WTA=ALPHA-ALPHA 1 
IF! FLAP. LT. -5.0) 
IFIFLAP.LT .-5.0 ) 
IF! FLAP. GT. 5.0) 

I FI FLAP .GT .5.0) 



7 

8 



10 

11 

12 

13 

14 
1 5 



16 



17 

18 
19 



20 



21 

22 



WRI TEI6, 100) 
ALFA=-6. 0 



ALFACtM 



GO 

GO 



TO 

TO 



FLAP 



FORMAT! «0' ,’ ERROR. FLAP 



WRITE!6» 101 ) 

FLAP=-5.0 
WRITEI6,101) FLAP 
FLAP= 5 . 0 

EXCEEDED LIMITS, 

IF! FLAP )4 , 5 , 6 
I B= 1 

WTB=- ( ( 5. O+FLAP ) / 5 . 0 ) 

GO TO 7 
I 8=2 
WTB=0. 0 
GO TO 7 

I 

WTB= I 5.0-FLAPJ/5.0 
IFIWTA) 8,9,10 
I A2= I A- 1 
GO TO 11 
I A2=I A 
GO TO 11 
I A2= I A+ 1 
I FI WTB ) 12,13,12 
I B2=2 
GO TO 14 
I B2= I B 

GO TCI 15, 16,1 
C0EF1=C0EFMX I 
C0EF2=C0EFMX I 
GO TO 18 

C0EF1 =COEFCX I IB, I A) 

C0EF2 = C0EFCX I 132, I A) 

GO TO 18 

C0EF1 =COEFOX I IB, I A) 

C0EF2=C0EF0X( IB2, I A) 

C0EF3= I COEF2-COEF1 )■ 

GO TO I 19,20,21 ) , M 
COEF 1 P=COEFMX I 13, IA2) 

C0EF2 P= COE FMX I I B2 , IA2) 

GO TO 22 

CCEF1P=C0EFCXI IB, IA2) 

C0EF2 P = COE FCX I IB2, IA2 ) 

GO TO 22 

COEF 1 P = COE FOX I I 3, I A2) 

COEF2P=COEFOX( 162, I A2 ) 

C0EF3 P= IC0EF2P-C0EF1P) *ABS l WTB ) +C0EF1P 
COEF=I C0EF3P-C0EF3) =i=ABS( WTA) +C0EF3 
RETURN 
END 

FORTRAN 



FLAP = ' 



7) ,M 
IB, IA) 
162, I A) 



: ABSI WTB ) +COEF 1 



33 ) 

, E 1 4 . 8 , 



, E 1 4 . 8 ) 
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SUBROUTINE INTRP1 



C 
C 
C 

c 

SUBROUT I NE INTRP1 1 BETA I , SOB t COEF S f M, THE T A , T HETA1 ) 
COMMON JUNK,CFDS1 ,CADS1 

DIMENSION JUNK! 297) ,CFDS1(1 1,11 )tCADSl( 11,11) 

DATA CONDEG/57. 295773/ 

SUB=SOB 

BETA=BETAl*CONDEG 

IFIBETA.GT.5.0) WRITE(6,100) BE T A 1 , M 
I F (BETA .GT .5. 0) BETA=5.0 

100 FORMAT! 'O' ,• ERROR. BETA EXCEEDS L I MI TS . BET A= • , E 14 . 8 , 2X, 
l'M=* , II ) 

I F(BETA .LT.-5.0) WRITE(6,100) BETA 1 » M 
IF(BETA.LT .-5.0) BETA=-5.0 

IF(BETA.LT.O.O) GO TO 1 
IF(BETA.GT.O.O) GO TO 2 
I A=6 
WTA=0 . 0 
GO TO 3 

1 I A1=-BETA 
I A=6- I A 1 
BET AG=- I A1 
WTA=BETA-BETAG 
GO TO 3 

2 I A1=BET A 
I A = 6 + I A 1 
BETAG=I A1 
WTA=BETA-BETAG 

3 IF ( SUB . LT .0.0) SU3=0.0 
IF(SUB.GT.IO.O) SUB=10.0 
IF! SUB . GT .0.0) GO TO 5 

I S= 1 
WTS=0 . 0 
G.O TO 7 
5 I S 1=SUB 
I S=IS1+1 
SU31=IS 1 
WTS= SU6-SUB1 

7 IF! WTA) 8,9,10 

8 I A2=I A- 1 
GO TO 11 

9 I A2=I A 
GO TO 11 

10 IA2=I A+ 1 

11 IF(WTS.GT.O.O) GO TO 12 
I S2= I S 

GO TO 15 

12 IS2=IS+1 

15 IF(M.EQ.2) GO TO 16 
C0EF1=C ADS1 ( IS , IA ) 

C0EF2 = CADS1( IS2, I A) 

GO TO 18 

16 C0EF1=CFDS1 ! IS, IA) 

C0EF2=CFDS1! IS2, I A) 

18 C0EF3=( C0EF2-C0EF1 )*WTS+C0EF1 
IF4M.EQ.2) GO TO 20 

19 COEF 1 P= CADS 1 ! IS, I A2) 

C0EF2P=CADS1 ! I S 2 , IA2) 

GO TO 22 

20 COEF 1P = CFDS 1 ! IS, IA2) 

C0EF2P=CFDS1! IS2, IA2) 

22 C0EF3P= (C0EF2P-C0EF1P) *WTS+CCEF1P 
C0EFS=(C0EF3P-C0EF3) ; '''ABS(WTA) +C0EF3 
THETA1=THETA*C0NDEG 
RETURN 
END 

//L.SYSLIB DO DISP=SHR,DSN=SYS1 .DSL , 

// VOL = S E R= MARY ,UNIT = 2314 

// DD D I SP = SHR ,DSN = SYS1 .FORTH B 

// DD D ISP=SHR,DSN=SYS1 .MPSLIB 
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//G.FT06F001 DD SP ACE = ( C YL , ( 3 ) ) 

//G.FT03F001 DO DSN = S 1 269 .GLOP , UN IT=2 32 1 ,DI SP = SHR , 

// VOL=SER=CEL004»LABEL=( » » * IN) 

//PLOT . SYSIN DD * 

STEP RESPONSE OF PITCH ANGLE 
NANAGARA 

0.0 2.0 - 2.0 1 . 5.0 6.0 
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